Stability of Leap-Frog Constant-Coefficients Semi-Implicit Schemes for the Fully Elastic System of Euler Equations. Case with Orography

The stability of onstantoe ients semi-impli it s hemes for the hydrostati primitive equations and the fully elasti Euler equations in presen e of expli itly treated thermal residuals has been theoreti ally examined in the earlier literature, but only for the ase of a at terrain. This paper extends these analyses to the ase where an orography is present, in the shape of a uniform slope. It is shown, with mass-based oordinates, that for the Euler equations, the presen e of a slope redu es furthermore the set of the prognosti variables whi h an be used in the verti al momentum equation to maintain the robustness of the s heme, ompared to the ase of a at terrain. The situation appears to be similar for systems ast in mass-based and height-based verti al oordinates. Still for mass-based verti al oordinates, an optimal prognosti variable is proposed, and shown to result in a robustness similar to the one observed for the hydrostati primitive equations system. The prognosti variables whi h lead to robust semi-impli it s hemes share the property of having umbersome evolution equations, and an alternative time-treatment of some terms is then required to keep the evolution equation reasonably simple. This treatment is shown not to modify substantially the stability of the time-s heme. This study nally indi ates that with a pertinent hoi e for the prognosti variables, mass-based verti al oordinates are equally suitable as height-based oordinates for e iently solving the ompressible Euler equations system. 2 1 Introdu tion The semi-impli it (SI) te hnique was introdu ed in meteorology by Robert et al. (1972), in order to in rease the numeri al e ien y with respe t to expli it s hemes, by allowing larger time-steps. SI s hemes are based on an arbitrary separation of the evolution terms between linear ontributions, treated impli itly, and non-linear (NL) residuals, treated expli itly. As dis ussed in Bénard et al., 2004 (BLSV04 hereafter), several lasses of SI s hemes an be de ned. In order to alleviate the problems linked to NL residuals, SI s hemes for whi h the linear terms have non horizontally-homogeneous oe ients (Thomas et al. 1998), or even nononstant oe ients (Skamaro k et al., 1997) an be designed. This latter type of s hemes is not examined here sin e the fo us of the present paper is ex lusively restri ted to the soalled lass of " onstantoe ients" SI s hemes, that is, to SI s hemes in whi h all the impli it linear terms have their oe ients onstant in time and horizontally homogeneous. For this lass of SI s hemes, the relatively large magnitude of NL residuals an result in instabilities, espe ially for long time-step. Note however that the stability an then be restored through an iteration of the s heme, as shown in Bénard, 2003 (B03 hereafter). A theoreti al framework for studying analyti ally the stability of these onstantoe ients SI s hemes in presen e of thermal NL residuals and in simpli ed ontexts has been proposed by Simmons et al. (1978), C té et al. (1983), and Simmons and Temperton (1997) for the hydrostati primitive equation (HPE) system, and extended to the fully elasti Euler equation (EE) system by B03. A more detailed review on the history of these stability analyses for SI s hemes an be found in this latter paper. Using the same theoreti al framework as in B03, a large sensitivity of the three-time levels (3-TL) 3 SI s heme stability to the hoi e of the prognosti variables was demonstrated for the EE system in BLVS04. In parti ular, the hoi e of the two new prognosti variables appearing in the EE system due to the relaxation of the hydrostati hypothesis was shown to have a dramati impa t on the stability. In the ase of a general mass-based oordinate (e.g. Laprise, 1992), an optimal set for these two variables was proposed as follows: P = p (1) d = p mRT w (2) where T is the temperature, w the verti al velo ity, p the true pressure, the hydrostati pressure, andm = ( = ) (in all this paper, see Appendix A for the de nition of symbols whi h are not de ned in the text). Finally, Bénard (2004) showed that an intrinsi instability of the two-time levels SI s heme for the EE system ould be eliminated by hoosing a slightly modi ed SI linear referen e system, whi h then an no longer be de ned as the tangent-linear operator of the omplete system around an existing referen e state. A ommon point to all the abovementionned studies was to negle t the orography as a sour e of nonlinear terms, and to fo us mainly on the instability indu ed by the thermal nonlinear residuals (as shown in BLVS04, the nonlinear terms indu ed by the pressure eld an be eliminated in mass-based oordinates by using appropriate oordinates, olumnintegrated mass variables, and pressure variables). Conversely, in the ase where all NL thermal residuals are negle ted, Ikawa (1988) showed that expli itly treated orographi terms ould make the evolution of ompressible systems in height-based oordinates unstable for the parti ular lass of prognosti variables and SI s hemes that he examined. It an therefore be suspe ted that a ombination of expli itly treated thermal residuals and 4 orographi terms ould result in new instabilities. In this ase, for mass-based oordinate systems, there is no proof that the hoi e (1) (2) proposed above is still the most relevant. This latter statement was on rmed experimentally, by showing that this hoi e a tually led to instability when an orography was introdu ed in the two-dimensional (verti al plane) version of the ooperative model "Aladin-NH", des ribed in Bubnová et al., 1995, (hereafter BHBG95). The instability ould be reprodu ed even for very simple ows, in luding initially balan ed resting isothermal ows over an isolated mountain. This experimental fa t prompted us to study more in detail the behaviour of the 3-TL SI EE system from the theoreti al point of view in presen e of orography, in order to investigate the nature of the asso iated instability and to seek possible remedies. The results of this study are reported in this paper. The stability analyses presented here are thus valid for the Aladin-NHmodel (BHBG95) as well as for numeri al models that would be based on the same prin iples (mass-based oordinates, similar linear SI separation, et ). In se tion 2, we will de ne an a ademi (simpli ed) framework whi h allows algebrai ally tra table stability analyses for the 3-TL SI EE spa eontinuous system in presen e of a simple orography whi h onsists in a "uniform slope mountain"; then the stability of the previously proposed hoi e (1) (2) will be examined in se tion 3. In se tion 4, an alternative variable dl will be proposed, and shown to result in a better stability in presen e of orography. In se tion 5, a numeri al assessment of the validity of the theoreti al results will be presented. In se tions 6 and 7 some omments on erning the HPE system and an alternative time-treatment of the verti al momentum equation will be developed, and the general on lusion will be presented in se tion 8. 5 2 Theoreti al framework for analyses The theoreti al framework used for the analyses in this paper is basi ally the same as in BLVS04, ex ept that an orography is introdu ed: the EE system is thus ast in the pure unstret hed terrain-following hydrostati -pressure-based = ( = s) oordinate and the ow is assumed to be dry, adiabati and invis id, in a non-rotating atmosphere and a Cartesian framework. The general set of equations (3)-(9) of BLVS04 is thus still valid, even in presen e of orography, and will be used here as a starting point. The general method for the analyses is also similar to B03 and BLVS04: an "a tual" steady basi -state X and a SI-referen e state X are hosen, and the analysis is performed for small perturbations around X (see B03 for deeper explanations and notations). Both X and X are assumed isothermal, resting and hydrostati ally-balan ed. The perturbations around X are assumed to remain small, hen e the sour e terms of the system an be linearized around X , and symboli ally noted L. Similarly, the sour e terms of the linear referen e system asso iated to the SI s heme are noted L . The SI s heme is then implemented a ording to Eq. (12) of BLVS04, and the residual (L L ) is treated expli itely. The domain is taken as a two-dimensional verti al plane along the x dire tion. The temperature of the basi a tual state T deviates from the one of the SI-referen e state T , thus generating expli itly treated thermal residuals. As outlined above, the major di eren e between the framework used in BLVS04 and the present one is that a uniform slope s is now assumed for the surfa e height zs of the bottom boundary: zs(x) = s x; (3) 6 hen e the surfa e geopotential is given by s(x) = gsx. Sin e the domain is assumed unbounded in the x dire tion in order to later allow simpler normal mode analyses, some are must be taken to assess the physi al relevan e of this framework. The olumn-integrated mass variations are des ribed through the prognosti variable q = ln( s). The stationarity of X for the omplete system imposes the pressure-gradient for e to vanish in the horizontal momentum equation. Hen e (using the above assumptions): RTrq +r s = 0 (4) or: rq = gs RT (5) where q(x) is the q eld in the a tual state X and r is the ( = x) operator along onstant surfa es. Under all mentionned assumptions, it an be he ked easily that X is then a stationary state for the omplete system (3) (9) of BLSV04, and onsequently for the linearized system asso iated to L as well. It appears that both elds s and q are unbounded when x be omes unbounded. However, this mathemati al feature does not indu e any parti ular problem from the physi al point of view, sin e only spatial or temporal variations of these elds appear in the equations. In other words, even if the height of the orography is unbounded along x, all orographi sour e terms always remain bounded sin e the slope s is a nite number. From the physi al point of view, the proposed framework is thus perfe tly relevant for des ribing the ow over a slanted unbounded orography. As outlined in B03, due to the elimination of upper and lower boundary onditions whi h will be performed in all these analyses, a more "lo al" point of view an be adopted 7 and the framework an then be onsidered as relevant also for des ribing the evolution of small-s ale perturbations inside a limited region of the atmosphere when the larger s ale environment is given by X . Using this lo al point of view, the slope s must therefore be understood as the mean slope of the orography at a s ale mu h larger than the s ale of the perturbations onsidered in the analysis. In this paper, we will examine atmospheri perturbations at the kilometri horizontal s ale, onsistently with future targets for NWP appli ations, hen e, the slopes retained in the analyses will be onsistent with terrestrial slopes at an horizontal s ale of 50-100 km, that is, s 2 [ 0:05;+0:05℄. For smaller (i.e. he tometri ) perturbation s ales, steeper slopes would have to be onsidered. The small perturbation f X around the X state is de ned by: e T = T T (6) e q = q q: (7) Sin e the other prognosti variables have a zero referen e-value, the tilde symbol is omitted for them. For the spe i ation of q in the X state, we follow the approa h usually adopted in pra ti al NWP appli ations, that is, q = ln( 00), where 00 is an arbitrary onstant. The omplete set of equations is then linearized around this q value, assuming no orography. 3 Stability analysis with the fP ; dg set of variables In this se tion the stability of the 3-TL SI EE system is examined for the set of variables fP; dg whi h was proposed in Eqs. (22) and (65) of BLVS04 for eliminating the problems linked with a dis repan y between T and T in the ase of at terrain. The linearization of the system is performed in the same way as in BLVS04, but now retaining all linear 8 terms involving s through r s or rq. The linearized expression of the three-dimensional divergen e D3 is: D3 = D + d + gs RT ~ u (8) where u is the horizontal wind omponent, D = ru, and ~ = ( = ). Additionally, the expression of ( _ = ) is required for the equation of P sin e: dP dt = (1 + P) _ pp _ ! (9) The linearized version of ( _ = ) writes: _ = SD gs RT (I S)u (10) and the linearized L system with orography be omes: D t = RTrr0 "G e T T (G I)P# RTr2e q (11) d t = g2 RT ~ (~ + I)P (12) e T t = RT Cv (r0u+ d) (13) P t = Sr0u Cp Cv (r0u+ d) (14) e q t = ND + gs RTNu (15) where the operator r0 is de ned by: r0 = r+ s H ~ (16) and H = (RT=g) (see Appendix A for other notations). The SI linear system L is obtained dire tly from the above system (11) (15) by substituting T ! T and s ! 0 9 (and onsequently r0 !r); it is of ourse not modi ed with respe t to the ase without orography. As a onsequen e, all terms asso iated to s are treated expli itly in (11) (15). The method of analysis then follows the one proposed in B03, and the reader is invited to refer to this paper for the details of notations and the algebrai developments. First, in this paragraph, the above system is shown to ful l the four onditions [C1℄ [C4℄ de ned in B03, and whi h are required for making possible the spa e ontinuous analyses with the proposed method. The number of prognosti variables is P = 4 in the sense of B03, and the spa e ontinuous state-ve tor is X = (X1; : : : ;X4) = (D; d; e T ;P). The linear operator l in [C1℄ involves l1 = ~ , l3 = r, and l4 = (~ + I)r, respe tively applied to (11), (13), and (14), as in se tion 7 of BLVS04. The "unbounded" linear system then be omes: ~ D t = RTrr0 " e TT (~ + I)P# (17) d t = g2 RT ~ (~ + I)P (18) r e T t = RT Cv (r0D +rd) (19) (~ + I)r P t = r0D Cp Cv (~ + I)(r0D +rd) (20) The stru ture equation, whi h allows to determine the timeontinuous normal modes of this linear system writes: 1 2 4 t4 + 2 t2 r02 + ~ (~ + I) H2 !+N2r02 = 0 (21) where: 2 = Cp CvRT (22) 10 N2 = g2 CpT (23) This stru ture equation appears to be formally similar to the stru ture equation without orography in BLVS04, ex ept that the original horizontal gradient operator r is repla ed by r0. The determination of the ontinuous normal modes thus follows dire tly from the ase without orography: the ondition [C'2℄ requires T > 0 and the stru ture of the normal modes of the ontinuous system is then given by: Xj(x; ) = Xj exp (ikr + s 2H )x (i 1=2) for j 2 (1; : : : ; 4) (24) where Xj is the omplex magnitude for the onsidered variable, and (kr; ) 2 IR. Note that is a non-dimensional verti al wave-number, and = 2 represents a mode with a verti al wavelength equal to the hara teristi height H of the atmosphere (here and in later similar uses related to waves geometry, is of ourse 3.1415...). The real value =kr represents the horizontal half-wavelength of the mode, that is, the distan e between two onse utive zeros of the real part of the omplex mode along the x dire tion. For su h a linear pertubation Xj(x; ), the energy density of the perturbation is proportional to X 2 j where = ( =RT ) (see e.g. Bannon, 1995). This energy density de omposes itself in three parts: kineti , available potential, and available elasti energy density. The normalmodes (24) have an horizontal variation along onstant surfa es whi h is onsistent with the growth of the modes with height due to the Boussinesq e e t and with the elevation of the terrain along the x axis. The normal modes Xj thus have a uniform amplitude along true horizontal surfa es only. However, the energy density of these normal modes is spatially uniform, as in the at-terrain ase, due to the ompensating variation of the mass density : 11 r( X 2 j ) = ~ ( X 2 j ) = 0 (25) For an eigenmode hara terized by (kr; ), the eigenvalues of derivative operators are: r = ikr + s 2H = ik (26) ~ = i 12 (27) r0 = i kr + s H = ik0: (28) The veri ation of [C3℄ [C4℄ pro eeds easily, as in the ase without orography. It should be noted that for the referen e system lL , ther operator is used everywhere, thus leading to the eigenvalue found in (26). For [C3℄, we have: 1 = (i 1=2) (29) 3 = ik (30) 4 = ik(i + 1=2); (31) and for [C4℄: 13 = Rkk0 (32) 13 = Rk2 (33) 14 = RTkk0(i + 1=2) (34) 14 = RT k2(i + 1=2) (35) 24 = g2 RT ( 2 + 1=4) (36) 24 = g2 RT ( 2 + 1=4) (37) 12 31 = RT Cv ik0 (38) 32 = RT Cv ik (39) 31 = 32 = RT Cv ik (40) 41 = ik0 1 Cp Cv (i + 1=2) (41) 41 = ik 1 Cp Cv (i + 1=2) (42) 42 = = 42 = ikCp Cv (i + 1=2): (43) All other ij oe ients are zero, and [C1℄ [C4℄ are nally ful lled. The analysis then pro eeds as for the at-terrain ase in BLVS04: For the 3-TL SI s heme, a numeri al growth-rate is introdu ed through: X (t) = X (t t) (44) X (t+ t) = 2X (t t); (45) (46) and the stability equation an be expressed as in BLVS04: Det(M) = 0; (47) where M is given by (46) (49) of BLVS04, used with the above values of j and ij. This eighth-degree stability polynomial equation in an be solved numeri ally: for any geometri al stru ture de ned by a pair (kr, ), the modulus of the eight roots 1(kr; ); : : : ; 8(kr; ) gives the growth rate of the eight orresponding eigenmodes. The growth-rate of the geometri al stru ture (kr, ) is then de ned by the maximum modulus of the eight roots: 13 (kr; ) = Max [ 1(kr; ); : : : ; 8(kr; )℄ (48) If one of the roots has a modulus larger than one, then the orresponding geometri al stru ture is unstable. In this paper, the asymptoti growth-rate for a given geometri al stru ture (kr, ) is de ned by the value of the above growth-rate in the limit of large time-steps. As dis ussed in B03, examination of asymptoti growth-rates is relevant sin e SI s hemes are used with long time-steps in pra ti e, and thus, asymptoti growth-rates provide a good indi ation of the robustness of a s heme independently of the parti ular value of the time-step. For onvenien e, a parameter for the thermal nonlinearity an be introdu ed through = T T T : (49) As an illustration of the results, the asymptoti growth-rate for the kilometri s ale horizontal mode with kr = 0:00157 m 1 and for T = 280 K is depi ted in Fig. 1 as a fun tion of and s. The growth-rate whi h is plotted is the maximum growth-rate obtained when repeating the above analysis for dis rete values of des ribing the interval [2 , 100℄. This interval represents verti al wavelenths varying between 500 m and H. In pra ti e however the maximum growth-rate for this gure is rea hed for the shortest verti al mode = 100 (not shown). The domain where the growth-rate is smaller than 1.1 is restri ted to a very small area along the axes, whi h means that for ows with signi ant values of the slope and the thermal nonlinearity, the s heme is highly unstable. 14 4 Stability analysis with the fP ; dlg set of variables As seen in BLVS04, the hoi e of the prognosti variables has a large impa t on the stability of the SI EE system, and the presen e of large nonlinear residuals in the elasti termD3 is suspe ted to potentially lead to instabilities. In presen e of orography, when d is used as prognosti variable, the elasti term D3 a tually has an expli itly treated thermal residual, given by the last RHS term of (8) in the linearized ontext. This residual, whi h is proportional to the slope, an thus be suspe ted to explain the instability found in the previous se tion when a slope is introdu ed. A new prognosti variable dl whi h avoids this problem is examined here. In the general mass-based hybrid terrain following oordinate of Laprise, 1992, this new variable dl is de ned by: dl = d + p mRT V :r (50) where is the geopotential. In pure oordinate, and using P variable, this be omes: dl = d + (1 + P) RT V :r (51) The elasti term D3 writes in all ases: D3 = D + dl (52) The new variable dl thus totally eliminates the thermal nonlinear residuals in the elasti term, even in presen e of slanted oordinate surfa es. The derivation of the linear system L with the variable dl is straightforward starting from (17)-(20). The geopotential gradient writes in oordinate: r = r s +R Z 1 r T 1 + P d 0 0 (53) 15 whi h in the urrent linear ontext yields: r = gs+RTGr e TT P! (54) The general relationship (50) between d and dl thus be omes in the linear ontext: dl = d+ gs RT ~ u (55) and the linear evolution equation for dl writes: dl t = d t + gs RT ~ u t (56) The original linearized system (11)-(14) for the d variable is thus modi ed into: D t = RTrr0 "G e TT (G I)P# RTr2e q (57) dl t = g2 RT ~ (~ + I)P + gsr0 " e TT (~ + I)P# (58) e T t = RT Cv (D + dl) (59) P t = Sr0u Cp Cv (D + dl) (60) e q t = ND + gs RTNu (61) Using the same operators l1, l3 and l4 as previously, the linear system is then modi ed into the following "unbounded" version: ~ D t = RTrr0 " e TT (~ + I)P# (62) dl t = g2 RT ~ (~ + 1)P + gsr0 " e TT (~ + I)P# (63) r e T t = RT Cv (rD +rdl) (64) 16 (~ + I)r P t = r0D Cp Cv (~ + I)(rD +rdl) (65) The stru ture equation is still given by (21), the stru ture of normal modes by (24), and the eigenvalues of spatial operators by (26) (28). The i oe ients are un hanged, and the ij oe ients whi h are modi ed with respe t to the previous analysis with d are: 23 = gs T ik0 (66) 24 = g2 RT ( 2 + 1=4) gsik0(i + 1=2) (67) 31 = RT Cv ik (68) 41 = ik0 Cp Cv ik(i + 1=2) (69) These new oe ients are used to build the new matrix M and the analysis pro eeds similarly to the previous ase, by solving (47). Figure 2 shows the asymptoti growth-rate in the same onditions as for Fig. 1 but for the new variable dl instead of the variable d (it is worth noti ing that the gure has no reason to be symmetri around the s = 0 axis sin e the stability may depend on the relative value of the oordinate slope and of the wave-surfa es slope). The domain for whi h the growth-rate is smaller than 1.1 is now mu h larger than in Fig. 1, learly indi ating an enhan ed robustness of the s heme with the new prognosti variable dl for this stru ture. In opposition with the d variable, the maximum growth-rate depi ted on the gure is not always obtained for the shortest verti al mode: for instan e the unstable areas lo ated near s = 0:05 and 0:2 are not due to the shortest verti al mode (not shown). However, even if using the new variable dl globally enhan es the stability, it does not guarantee a formal stability of the s heme for arbitrarily long time-steps in the ontext of the analysis. As indi ated by the shape of the ontour line = 1:01 in Fig. 2, the s heme 17 is in fa t formally stable only for s = 0 and 2 [ 0:5; 1℄. As soon as s 6= 0, the dl variable still leads to a formal instability (although mu h weaker than for the d variable). This behaviour ontrasts with the formal stability obtained with d for at terrains in BLVS04, even in the long time-steps limit. The analyses presented above thus indi ate that when orography is present, the EE system is not likely to be solved in a stable way with a SI s heme using arbitrarily long time-steps even with the dl variable. The stability of the same horizontal mode kr = 0:00157 as in the previous gures is now examined when the time-step is redu ed to t = 100 s. This time-step represents a realisti (although ambitious) value for a model whi h would resolve the onsidered mode kr with a 3-TL SI time-dis retization. Figure 3 shows the growth-rate obtained for this time-step with the d variable (note the modi ed shading limits ompared to Figs. 1 and 2). Not surprisingly, the growth-rate de reases when the time-step is redu ed, however, the values of the growth-rate when s and signi antly deviate from 0 are still in ompatible with a pra ti al use. The growth-rate obtained in the same onditions with the variable dl is depi ted in Fig. 4. Here also, a formal stability is obtained only for s = 0 as indi ated by the shape of the ontour line = 1:001, however, the value of the growth-rate is very lose to 1 in a wide area (in pra ti e has a typi al absolute magnitude smaller than 0.2), and the s heme an be onsidered as viable when used in a spa e-dis retized model, in whi h many pro esses are likely to a t in a di usive fashion. 5 Numeri al assessment in the a ademi ontext In order to demonstrate the qualitative relevan e of the previous analyses, a similar a ademi ontext is reprodu ed in the spe tral Aladin-NH numeri al model do umented in BHBG95. A non-rotating adiabati verti al-plane version of the model is used, with a 18 resting, isothermal and hydrostati ally-balan ed initial state. The horizontal domain is horizontally y li , 64 km wide with a 21 spe tral trun ation, thus allowing a minimum horizontal wavelength of about 3048 m. The verti al domain onsists in 60 regularly spa ed levels. The value of the time-step is t = 50 s. For pra ti al reasons, the orography does not onsist of a onstant slope as in the analyses but is taken as a osine fun tion of same wavelength as the horizontal width of the domain. The deviation of the osine orography is 500 m, thus resulting in two areas of opposed maximum slope s = 0:05. The model is made free of any damping pro ess (horizontal di usion, SI timelter,...), however, a weak lassi al Asselin timelter ( = 0:01) is applied to prevent a separation of physi al and omputational modes resulting from the leap-frog time-dis retization. A very small random initial wind perturbation is introdu ed in order to prevent an undeterministi evolution when the initial equilibrium state is physi ally stable but numeri ally unstable. The atmospheri and referen e temperature are set to 285 K and 350 K respe tively (hen e -0.186). The experiment with d numeri ally diverges after 49 time-steps: an unstable ow develops above the areas of maximum slope. For the experiments with dl, the integration remains stable after 2000 time-steps: the dis retization pro esses a ting in a numeri al model (in luding semi-Lagrangian interpolations) are able to stabilize unstable modes when their growth-rate is very small. However, the instability for large time-steps appears when t is pushed to 200 s: the experiment with d and dl respe tively diverges after 10 and 41 time-steps, qualitatively on rming the results of the previous analyses. 19 6 Comments Sin e no stability analysis of onstantoe ients SI s hemes with orography and NL thermal residuals has been reported in the earlier literature, it is interesting to examine what is the situation for the HPE system. A similar (although more simple) analysis an be performed for the HPE system. In this ase, the L system redu es to: D t = Rrr0G e T RTr2e q (70) e T t = RT Cp Sr0u (71) e q t = ND + gs RTNu: (72) Following the formalism of B03, a linear verti al operator hl1 = ~ ; l2 = (~ + I)ri an be applied to the previous system: the number of prognosti variables for the unbounded system is then P = 2, the state variable is (D, e T ). For the [C3℄ ondition, we have: 1 = (i 1=2) (73) 2 = ik(i + 1=2); (74) and for [C4℄: 12 = Rkk0 (75) 12 = Rk2 (76) 21 = RT Cp ik0 (77) 21 = RT Cp ik (78) (79) 20 All other ij oe ients are zero. The analysis then follows as in previous se tions. The growth-rate for t = 100 s and kr= 0.00157 as in Figs. 3 and 4 is depi ted in Fig. 5 for the HPE system. A weak instability is present in the onsidered domain with a magnitude roughly similar to the one obtained with the dl variable. Hen e the weak instability of the SI s heme in presen e of orography found with dl is not a new feature of the EE system, but was already present in the HPE system, although not reported in the earlier literature. This may be explained by the fa t that this instability is maybe too weak to signi antly endanger urrent pra ti al NWP appli ations. It should be noted that the dramati sensitivity to the hoi e of the prognosti variables with sloped-terrain, as dis ussed in se tions 3, 4, also exists for EE systems in heightbased oordinates with onstantoe ients SI s hemes, su h as e.g. in the CRCM model (Caya and Laprise, 1999, CL99 hereafter): as dis ussed in BLVS04, the natural hoi e for the verti al momentum variable for these models would be the verti al velo ity w, but this hoi e was experimentally found to result in a very unstable SI s heme with sloped terrain (Laprise, personal ommuni ation). To solve this problem, these models use, as a prognosti variable, the ontravariant verti al velo ityW [given by (21) of CL99℄, whi h, in the ontext of height-based oordinated is the ounterpart of dl in mass-based oordinates. 7 Alternative time treatment for the ross-term of dl The use of dl in mass-based oordinates, as well as W in height-based oordinates (when a onstantoe ients approa h is used as in CL99), leads to a pra ti al problem for the omplete (nonlinear) system, sin e these prognosti variables have rather umbersome evolution equations. Considering (50), the evolution equation for dl would involve a very 21 large number of terms oming from the evolution of the soalled ross-term, that is, the last term in the right hand side (rhs) of (50). Similar ompli ations o ur for W in height-based Gal-Chen oordinates. The solution adopted to ir umvent this problem is illustrated here for the dl variable in the general hybrid mass-based oordinate , but a similar argument an be developped for W in height based oordinates. The dl equation is written as follows: ddl dt = RHS dd dt!+ d dt p mRT V :r ! ; (80) where RHS (dd=dt) is the rhs of the omplete (nonlinear) evolution equation for d. This latter rhs is mu h simpler than the total rhs of (80), and is treated in the lassi al SI fashion: it is separated into a linear part L d(X ) and a non-linear residualMd(X ) L d(X ) exa tly as for the d equation. Sin e the oe ients ij are identi al in se tions 3 and 4, the linear part L dl(X ) of the dl evolution equation is equal to L d(X ). Then the evolution equation for dl is time-dis retized as: dl+ dl 2 t = Md(X 0) L d(X 0) + L d(X+) + L d(X ) 2 + 1 t 24 p mRT V :r !0 p mRT V :r ! 35 (81) where supers ripts , +, and 0 represent variables at time (t t), (t + t) and t respe tively. Hen e, the linear part of the dl equation has a entred impli it treatment, the (Md L d) part is treated expli itely, and the last term of (80) is treated in a diagnosti way from the two most re ent available states. The same approa h is applied in models with height-based verti al oordinates and onstantoe ients SI s hemes, in order to avoid the umbersome evolution equation of W . For instan e in CL99, the ombination of their Eqs. (48) and (52) leads to a single 22 prognosti equation for W formally similar to the above equation (81). The stability analysis for the SI s heme with the dl variable presented in se tion 4 an be easily adapted to take into a ount the spe i time-treatment of the dl equation. A l2 = r operator is applied to (56) in order to eliminate u in favour of D in the last term. Hen e for the [C3℄ ondition: 2 = ik (82) and two ij oe ients must be modi ed with respe t to se tion 4: 23 = 0 (83) 24 = g2 RT ik( 2 + 1=4): (84) A new set of oe ients ij must be introdu ed: 21 = gs RT (i 1=2); (85) all other ij oe ients being equal to zero. Finally, the matrixM1 in Eq. (47) of BLVS04 must be modi ed to: (M1)ij = Æij t ij i 2( 1) ij i (86) For t = 100 s, the results for the stability are not signi antly modi ed with respe t to those obtained with the normal treatment of dl as seen in Fig. 6. For long timesteps, this time treatment brings a further in rease of stability ompared to the normal treatment (not shown). These spe i time treatments with dl andW variables thus ombine the advantages of an in reased simpli ity and a similar robustness. Their main disadvantage is that they are 23 only rst-order a urate in time. The generally a epted requirement that NWP models should be globally se ond-order a urate in time is thus violated for the terms treated in this way. However, the impa t of these treatments an be easily evaluated in a given model by hanging the size of the time-step, and no parti ular problem linked to this point has been reported so far (e.g. in the CRCM model). Moreover, it should be noted that this disdvantage disappears by onstru tion for iterative entered impli it (ICI) s hemes as de ned in B03. These ICI s hemes, based on a larger number of iterations of the linear impli it equation inversion, have been shown in B03 to be more robust than the SI s heme, and ould be an interesting alternative for the fully ompressible meso-s ale NWP in the future, espe ially when the the physi al parameterizations are ostly ompared to the dynami al kernel of the onsidered model. 8 Con lusion The stability of onstantoe ients 3-TL SI s hemes for the EE and HPE systems has been examined in the ase of a uniform slope orography in presen e of expli it thermal residuals. It has been shown that for the EE system, the presen e of this slope redu es furthermore the set of the prognosti variables whi h an be used for the verti al momentum equation. The only possible variables are those for whi h the elasti term D3 has no nonlinear residual in presen e of slope. The variable proposed in the massoordinate ontext of this paper obviously meets this riterion, sin e it writes dl = D3 D, and D is also a prognosti variable. The resulting robustness has been shown to be similar to the one observed for the HPE system in presen e of slope. For prognosti variables su h as dl (in mass-based oordinates) or W (in height-based oordinates), the evolution equation is umbersome, and an alternative time-treatment of some terms is required. However, 24 further analyses have shown that this treatment does not modify signi antly the stability of the SI s heme. This paper, together with B03, BLVS04 and B04 forms a omplete set whi h allows to better understand the behaviour of the EE system with onstantoe ients SI (or ICI) time-s hemes from the theoreti al point of view. The general learning drawn from this set of papers is that provided a great are is taken about all the details of the formulation, the robustness of the EE system with onstantoe ients SI s hemes should be omparable to the one of the HPE system, and ompatible with meso-s ale NWP purposes. The modi ations proposed in this series of papers have been implemented in the Aladin-NH model (in luding the possibility of using ICI s hemes), resulting in what may nally be viewed as a new dynami al kernel for this model. Another general learning is that there is a deep formal dualism between the systems ast in height-based and mass-based oordinates as far as onstantoe ient SI s hemes are on erned. Ex ept the fa t that the spe i ation of a referen e pressure pro le is mandatory in height-based systems with onstantoe ients SI s hemes, a feature whi h has no equivalent in mass-based systems, no signi ant di eren e was found between the two lasses of systems in terms of behaviour and stability. Another illustration of this dualism an be found in the very simple form that takes the relaxation of the shallowatmosphere approximation in the mass-based EE system (Wood and Staniforth, 2003) similarly to height-based systems. We laim that this dualism invalidates a ommon (though not published) belief a ording to whi h height-based oordinates are the best suited ones to build e ients nonhydrostati models. Finally, it should be noted that the onstantoe ients s hemes examined here, when oupled with iterative (ICI) s hemes represent an alternative approa h to s hemes with nononstant oe ients as proposed by Skamaro k et al. (1997) or Thomas et al. (1998). 25 Together with C té et al. (1998) and Cullen (2001), we believe this approa h is worth being onsidered for solving the EE system in a stable and a urate manner. A knowledgments: Part of the resear h reported in this paper was supported by the ALATNET grant HPRN-CT-1999-00057 of the European Union TMR/IHP Programme. We gratefully a knowledge helpful dis ussions with Jozef Vivoda. 26 Appendix A : List of Symbols verti al spatial operators in oordinate: GX = R 1 (X= 0)d 0 SX = (1= ) R 0 Xd 0 NX = R 1 0 Xd IX = X ~ X = ( X= ) mis ellaneous symbols: r : ( = x) along onstant levels of the onsidered verti al oordinate. r : ve tor horizontal gradient operator along onstant levels of the onsidered verti al oordinate. ( = t) : Eulerian time derivative (d=dt) : Lagrangian time-derivative g : gravitational a eleration R, Cp, Cv : dry air thermodynami onstants s : hydrostati surfa e-pressure (proportional to the olumn-integrated mass in a Cartesian system as outlined in Wood and Staniforth, 2003). 27 Referen es Bannon, P. R., 1995: Hydrostati adjustment; Lamb's problem. J. Atmos. S i, 52, 1743 1752. Bénard, P., 2003: Stability of semi-impli it and iterative entered-impli it time disretizations for various equation systems used in NWP. Mon. Wea. Rev., 131, 2479-2491. Bénard, P., 2004: On the use of a wider lass of linear systems for the design of onstantoe ients semi-impli it time-s hemes in NWP. Mon. Wea. Rev., 132, 1319-1324. Bénard, P., R. Laprise, J. Vivoda and P. Smolíková, 2004: Stability of the leap-frog onstantoe ients semi-impli it s heme for the fully elasti system of Euler equations. Flat terrain ase. Mon. Wea. Rev., 132, 1306-1318. Bubnová, R., G. Hello, P. Bénard, and J.F. Geleyn, 1995: Integration of the fully elasti equations ast in the hydrostati pressure terrain-following oordinate in the framework of the ARPEGE/Aladin NWP system. Mon. Wea. Rev., 123, 515-535. Caya, D., and R. Laprise, 1999: A semi-impli it semi-Lagrangian regional limate model: the Canadian RCM. Mon. Wea. Rev., 127, 341-362. C té, J., M. Béland, and A. Staniforth, 1983: Stability of verti al dis retization s hemes for semi-impli it primitive equation models: theory and appli ation. Mon. Wea. Rev., 111, 1189-1207. C té, J., S. Gravel, A. Méthot, A. Patoine, M. Ro h, and A. Staniforth, 1998: The Operational CMC-MRB Global Environmental Multis ale (GEM) Model. Part I: Design Considerations and Formulation. Mon. Wea. Rev., 126, 1373-1395. 28 Ikawa, M., 1988: Comparison of some s hemes for nonhydrostati models with orography. J. Meteorol. So . Jpn., 66, 753-776. Laprise, R., 1992: The Euler equations of motion with hydrostati pressure as an independent variable. Mon. Wea. Rev., 120, 197-207. Qian, J.-H., F. H. M. Semazzi, and J. S. S roggs, 1998: A global nonhydrostati semiLagrangian atmospheri model with orography. Mon. Wea. Rev., 126, 747-771. Robert, A. J., J. Henderson, and C. Turnbull, 1972: An impli it time integration s heme for baro lini models of the atmosphere . Mon. Wea. Rev., 100, 329-335. Simmons, A. J., C. Temperton, 1997: Stability of a two-time-level semi-impli it integration s heme for gravity wave motion. Mon. Wea. Rev., 125, 600-615. Simmons, A. J., B. Hoskins, and D. Burridge, 1978: Stability of the semi-impli it method of time integration. Mon. Wea. Rev., 106, 405-412. Skamaro k, W. C., P. K. Smolarkiewi z, and J. B. Klemp, 1997: Pre onditioned onjugateresidual solvers for Helmholtz equations in nonhydrostati models. Mon. Wea. Rev., 125, 587-599. Wood, N. and A. Staniforth, 2003: The deep-atmosphere Euler equations with a massbased verti al oordinate. Quart. J. Roy. Meteor. So ., 129, 1289-1300. Thomas, S.J., C. Girard, R. Benoit, M. Desgagné, and P. Pellerin, 1998: A new adiabati kernel for the MC2 model. Atmos. O ean, 36 (3), 241-270. 29 List of Figures Fig. 1: Growth-rate for the simpli ed problem as a fun tion of the nonlinearity parameter and slope s for the variable d with t = +1. Fig. 2: Growth-rate for the simpli ed problem as a fun tion of the nonlinearity parameter and slope s for the variable dl with t = +1. Fig. 3: Growth-rate for the simpli ed problem as a fun tion of the nonlinearity parameter and slope s for the variable d with t = 100 s. Fig. 4: Growth-rate for the simpli ed problem as a fun tion of the nonlinearity parameter and slope s for the variable dl with t = 100 s. Fig. 5: Growth-rate for the simpli ed problem as a fun tion of the nonlinearity parameter and slope s for the HPE system with t = 100 s. Fig. 6: Growth-rate for the simpli ed problem as a fun tion of the nonlinearity parameter and slope s for the dl variable with the alternative time-treatment, and with t = 100 s. 30